home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / a_utils / expanded.lha / expanded / src / SubSet.C < prev    next >
C/C++ Source or Header  |  1992-03-19  |  46KB  |  1,592 lines
  1. //
  2. // Linear-Affine-Projective Geometry Package
  3. //
  4. // SubSet.C
  5. //
  6. // $Header$
  7. //
  8. // William J.R. Longabaugh 
  9. // University of Washington
  10. //
  11. // Implementation of the linear-affine-projective geometry
  12. // package described in William J.R. Longabaugh, "An Expanded
  13. // System for Coordinate-Free Geometric Programming", Master's 
  14. // thesis, University of Washington, 1992.
  15. //
  16. // Copyright (c) 1992, William J.R. Longabaugh
  17. //   Copying, use and development for non-commercial purposes permitted.
  18. //   All rights for commercial use reserved.
  19. //   This software is unsupported and without warranty; it is
  20. //   provided "as is".
  21. //
  22. // ***********************************************************************
  23.  
  24. #include <math.h>
  25. #include "Lap.h"
  26.  
  27. // ***********************************************************************
  28. //
  29. // Class for exclusive use of SubSet class to hold onto data for a 
  30. // subset.  What is really wanted is a simple structure; I am not
  31. // really interested in data hiding here.  However, under g++ 1.37.1,
  32. // this apparently MUST be a class, not a struct, so that virtual 
  33. // ptrs of class members are initialized.  Also, the fact that class
  34. // members will be automatically initialized using null constructors
  35. // is desirable.  Thus, this is a class, but note that all members
  36. // are public, allowing unlimited access.
  37. // Note that since there is a VSubSet in this block, SubSets
  38. // must in fact be implemented as pointers to information blocks:
  39. //
  40.  
  41. class SubSetInfo {
  42.  
  43.   public:
  44.  
  45.     SubSetType t;                // Type of subset
  46.     char     name[MAX_NAMESIZE];        // Name of subset
  47.     Space    s;                         // Space containing subset
  48.     int      d;                // Dimension of subset
  49.     Matrix   tstmtx;            // Testing matrix
  50.     Matrix   fullbas;            // Matrix rep. of full basis
  51.     int      isOffset;            // Column for testing
  52.     int      substart;            // Column for testing
  53.     int      zerostart;            // Column for testing
  54.     Scalar   offset;            // Offset of affine hyperplane
  55.     Matrix   fromfull;            // Matrix from full to subset basis
  56.     Matrix   tofull;            // Matrix from subset to full basis
  57.     GeObType accepts;            // Type of object in subset
  58.     VSubSet  tansub;            // Tangent subset
  59.  
  60. // Constructors/Destructors:
  61.  
  62.         SubSetInfo(void)  {}         // Null constructor
  63.         ~SubSetInfo(void) {}         // Destructor
  64.         SubSetInfo(char* namein, int d); // Used to build error info block
  65. };
  66.  
  67. // ***********************************************************************
  68. //
  69. // Declare the existence of the error info block:
  70.  
  71. static SubSetInfo errinfo;
  72.  
  73.  
  74. // ***********************************************************************
  75. // ***********************************************************************
  76. //
  77. // SubSetInfo Class
  78. //
  79. // ***********************************************************************
  80. // ***********************************************************************
  81. //
  82. // Special constructor used to build an "error info block".  Uninitialized
  83. // subsets point to that block.
  84.  
  85. SubSetInfo::SubSetInfo(char* namein, int dim)
  86. {
  87.   t = NULL_SUBSET;
  88.   d = dim;
  89.   isOffset = -1;
  90.   substart = -1;
  91.   zerostart = -1;
  92.   offset = 0.0;
  93.   accepts = NULL_GEOB;
  94.  
  95. // Local copy of the debug name:
  96.  
  97.   strncpy(name, namein, MAX_NAMESIZE - 1);
  98.   name[MAX_NAMESIZE - 1] = '\0';
  99. }
  100.  
  101. // ***********************************************************************
  102. // ***********************************************************************
  103. //
  104. // SubSet Class
  105. //
  106. // ***********************************************************************
  107. // ***********************************************************************
  108. //
  109. // Private/protected member functions
  110. //
  111. // ***********************************************************************
  112. //
  113. // Great candidates for inlining, but see the comment in the file
  114. // Space.C:
  115. //
  116. // ***********************************************************************
  117.  
  118. int SubSet::IsOffset(void) {return info->isOffset;}
  119.  
  120. // ***********************************************************************
  121.  
  122. Scalar SubSet::Offset(void) {return info->offset;}
  123.  
  124. // ***********************************************************************
  125.  
  126. Matrix& SubSet::FromFull(void) {return info->fromfull;}
  127.  
  128. // ***********************************************************************
  129.  
  130. Matrix& SubSet::ToFull(void) {return info->tofull;}
  131.  
  132. // ***********************************************************************
  133.  
  134. Matrix& SubSet::AugBasis(void) {return info->fullbas;}
  135.  
  136. // ***********************************************************************
  137. //
  138. // Used to normalize the representation of projective points contained
  139. // in an affine subset.  It does not do ANY error checks, so beware!
  140. //
  141.  
  142. GeOb SubSet::Normalize(GeOb& ob)
  143. {
  144.   Matrix newmat = ob.MatrixRep();
  145.   Scalar factor = info->offset / (newmat * info->tstmtx)[0][info->isOffset];
  146.  
  147.   GeOb retval(ob.Holds(), ob.SpaceOf(), factor * newmat);
  148.   return retval;
  149. }
  150.  
  151. // ***********************************************************************
  152. //
  153. // Given two affine subsets in projective spaces, determine the
  154. // scalar factor that must be applied to place the origin point
  155. // of the first at the offset of the secong origin point
  156. // It does not do ANY error checks, so beware!
  157. //
  158.  
  159. Scalar SubSet::NormVal(SubSet& targ)
  160. {
  161.   Matrix newmat = targ.ToFull()[targ.info->isOffset];
  162.   Scalar factor = info->offset / (newmat * info->tstmtx)[0][info->isOffset];
  163.   return factor;
  164. }
  165.  
  166. // ***********************************************************************
  167. // 
  168. // This function copies only a maximum number of characters into
  169. // the debug name buffer.
  170. //
  171.  
  172. void SubSet::CopyDebugName(char* buf)
  173. {
  174.   strncpy(info->name, buf, MAX_NAMESIZE - 1);
  175.   info->name[MAX_NAMESIZE - 1] = '\0';
  176.   return;
  177. }
  178.  
  179. // ***********************************************************************
  180. // 
  181. // This function builds a tagged debug name.  The user is responsible
  182. // for calling DestroyTagName() when finished.
  183. //
  184.  
  185. char* SubSet::BuildTagName(char* tag, Space& s) 
  186. {
  187.   int len = strlen(tag);
  188.   char* buf = new char[MAX_NAMESIZE + len];
  189.  
  190.   strncpy(buf, tag, len);
  191.   s.Name(buf + len);
  192.   return (buf);
  193. }
  194.  
  195. // ***********************************************************************
  196. // 
  197. // This function builds a tagged debug name.  The user is responsible
  198. // for calling DestroyTagName() when finished.
  199. //
  200.  
  201. char* SubSet::BuildTagName(char* tag, SubSet& s) 
  202. {
  203.   int len = strlen(tag);
  204.   char* buf = new char[MAX_NAMESIZE + len];
  205.  
  206.   strncpy(buf, tag, len);
  207.   s.Name(buf + len);
  208.   return (buf);
  209. }
  210.  
  211. // ***********************************************************************
  212. // 
  213. // This function kills a tagged debug name.
  214. //
  215.  
  216. void SubSet::DestroyTagName(char* buf) 
  217. {
  218.   delete buf;
  219. }
  220.  
  221. // ***********************************************************************
  222. //
  223. // This little routine is used by subset constructors to convert lists
  224. // of arbitrary geometric objects into objects of the desired type
  225. // in the desired space.  A matrix containing std. coords of the
  226. // objects is built:
  227. //
  228.  
  229. Boolean SubSet::ConvertList(GeObList& v, GeObType targ,
  230.                 Space& destspace, Matrix* temp)
  231. {
  232.   for (int i = 0; i < v.Length(); i++) {
  233.     GeOb hold = v[i].MapTo(targ);
  234.     if (hold.SpaceOf() != destspace) {
  235.       return (FALSE);
  236.     } 
  237.     (*temp)[i] = hold.MatrixRep()[0];
  238.   }
  239.   return (TRUE);
  240. }
  241.  
  242. // ***********************************************************************
  243. //
  244. // Modified Gram-Schmidt orthogonalization.  The "have" matrix has
  245. // orthonormal rows, and represents the existing orthonormal basis
  246. // of some subspace.  The "try" matrix contains vectors to use
  247. // for constructing a total of "dim" orthonormal vectors; any
  248. // vectors in the try matrix that are not linearly independent
  249. // are skipped.  We set the error flag if there are not enough vectors
  250. // to use for meeting the "dim" target:
  251. //
  252.  
  253. Matrix SubSet::GSO(Matrix& have, Matrix& try, int start,
  254.            int dim, Boolean* error)
  255. {
  256.   Matrix retval(dim, have.Columns());
  257.   int count = 0;
  258.  
  259. // Start by filling in with the vectors we already have:
  260.  
  261.   for (int i = 0; i < have.Rows(); i++) {
  262.     retval[i] = have[i];
  263.     count++;
  264.   }
  265.  
  266. // Go through the test vectors until we have enough:
  267.  
  268.   int tryrow = start;
  269.   while (count < dim) {
  270.     Matrix nextrow(1, have.Columns());
  271.     if (tryrow == try.Rows()) {
  272.       *error = TRUE;
  273.       return retval;
  274.     }
  275.     if (this->Orth(&retval, count, try, tryrow)) {
  276.       count++;
  277.     }
  278.     tryrow++;
  279.   }
  280.   *error = FALSE;
  281.   return retval;
  282. }
  283.  
  284. // ***********************************************************************
  285. //
  286. // Support for GSO.  Does the actual GSO for a vector.  Returns TRUE
  287. // iff the vector could be orthogonalized.  FALSE means it was
  288. // contained in the subspace we have built.  
  289. //
  290.  
  291. Boolean SubSet::Orth(Matrix* have, int count, Matrix& try, int tryrow)
  292. {
  293.   Scalar factor;
  294.   Scalar mag;
  295.  
  296. // Modified GSO does orthogonalization incrementally, once for each
  297. // vector we already have.  If we ever end up with the zero vector,
  298. // the attempt has failed and we need to return.  Otherwise, normalize
  299. // as we go:
  300.  
  301.   Matrix hold = try[tryrow];
  302.   for (int i = 0; i < count; i++) {
  303.     Matrix proj = (*have)[i];
  304.     factor = (proj * Transpose(hold))[0][0] / (proj * Transpose(proj))[0][0];
  305.     hold -= (factor * proj);
  306.     mag = sqrt((hold * Transpose(hold))[0][0]);    
  307.     if (fabs(mag) < EPSILON) {
  308.       return (FALSE);
  309.     } else {
  310.       hold = hold / mag;
  311.     }
  312.   }
  313.   (*have)[count] = hold[0];  
  314.   return (TRUE);
  315. }
  316.  
  317. // ***********************************************************************
  318. //
  319. // Build a subset from an info block:
  320. //
  321.  
  322. SubSet::SubSet(SubSetInfo* ts)
  323. {
  324.   info = ts;
  325.   type = ANY_SUBSET;
  326.   holds = ts->t;
  327. }
  328.  
  329. // ***********************************************************************
  330. //
  331. // Dumps data for a subset:
  332. //
  333.  
  334. void SubSet::data_out(ostream &c, int indent, char* name) 
  335. {
  336.   char *ibloc = new char[indent + 1];
  337.   for (int i = 0; i < indent; i++) {
  338.     *(ibloc + i) = ' ';
  339.   }
  340.   *(ibloc + indent) = '\0';
  341.  
  342.   c << ibloc << ast;
  343.   c << ibloc << name;
  344.   c << ibloc << "Type of subset: ";
  345.   SubSetTypeOut(c, type);
  346.   c << "\n";
  347.   c << ibloc << "Currently holds: ";
  348.   SubSetTypeOut(c, holds);
  349.   c << "\n";
  350.   if (holds != NULL_SUBSET) {
  351.     c << ibloc << "Name of subset: " << info->name << "\n";
  352.     c << ibloc << "Embedding space:\n";
  353.     info->s.debug_out(c, indent + ERR_IND);
  354.     c << ibloc << "Dimension: " << info->d << "\n";
  355.     c << ibloc << "Type of object in subset: ";
  356.     GeObTypeOut(c, info->accepts);
  357.     c << "\n";
  358.     c << ibloc << "Offset value column: " << info->isOffset << "\n";
  359.     c << ibloc << "Subset value starting column: " << info->substart << "\n";
  360.     c << ibloc << "Zero value starting column: " << info->zerostart << "\n";
  361.     c << ibloc << "Offset of affine hyperplane: " << info->offset << "\n";
  362.     c << ibloc << "Memory location: " << (long)info << "\n";
  363.     c << ibloc << "Matrix representation of testing matrix:\n";
  364.     info->tstmtx.debug_out(c, indent + ERR_IND);
  365.     c << ibloc << "Matrix representation of full basis matrix:\n";
  366.     info->fullbas.debug_out(c, indent + ERR_IND);
  367.     c << ibloc << "Matrix to convert to subset basis from full basis:\n";
  368.     info->fromfull.debug_out(c, indent + ERR_IND);
  369.     c << ibloc << "Matrix to convert to full basis from subset basis:\n";
  370.     info->tofull.debug_out(c, indent + ERR_IND);
  371.     c << ibloc << "Tangent subset:\n";
  372.     info->tansub.debug_out(c, indent + ERR_IND);
  373.   }  
  374.   c << ibloc << ast;
  375.  
  376.   delete ibloc;
  377.   return;
  378. }
  379.  
  380. // ***********************************************************************
  381. //
  382. // Public member functions
  383. //
  384. // ***********************************************************************
  385.  
  386. SubSet::SubSet(void) {type = ANY_SUBSET; holds = NULL_SUBSET; info = &errinfo;}
  387.  
  388. // ***********************************************************************
  389. //
  390. // Memberwise initialization:
  391. //
  392.  
  393. SubSet::SubSet(SubSet &s)
  394. {
  395.   info = s.info;
  396.   type = ANY_SUBSET;
  397.   holds = s.holds;
  398. }  
  399.  
  400. // ***********************************************************************
  401. //
  402. // Memberwise assignment:
  403. //
  404.  
  405. SubSet& SubSet::operator=(SubSet &s)
  406. {
  407. //
  408. // This weird checking is required to bypass the apparent inheritance of
  409. // assignment under g++ 1.37:
  410. //
  411.   if ((type != ANY_SUBSET) &&
  412.       ((type != s.holds) && (s.holds != NULL_SUBSET))) {
  413.     errh.ErrorExit("SubSet& SubSet::operator=(SubSet &)",
  414.                    "Illegal assignment attempted", *this, s);
  415.   }
  416.   holds = s.holds;
  417.   info = s.info;
  418.   return (*this);
  419. }  
  420.  
  421. // ***********************************************************************
  422. //
  423. // Output for debugging:
  424. //
  425.  
  426. void SubSet::debug_out(ostream &c, int indent) 
  427. {
  428.   static char* name = "SubSet object\n";
  429.   this->data_out(c, indent, name);
  430.   return;
  431. }
  432.  
  433. // ***********************************************************************
  434. //
  435. // Great candidates for inlining, but see the comment in the file
  436. // Space.C:
  437. //
  438. // ***********************************************************************
  439.  
  440. char* SubSet::Name(char *buf) {return (strcpy(buf, info->name));}
  441.  
  442. // ***********************************************************************
  443.  
  444. Space SubSet::EmbeddingSpace(void) {return (info->s);}
  445.  
  446. // ***********************************************************************
  447.  
  448. GeObType SubSet::Accepts(void) {return (info->accepts);}
  449.  
  450. // ***********************************************************************
  451.  
  452. int SubSet::Dim(void) {return (info->d);}
  453.  
  454. // ***********************************************************************
  455. //
  456. // Returns TRUE iff the subset spans the entire embedding space.
  457. //
  458.  
  459. Boolean SubSet::IsFullSpace(void)
  460. {
  461.   if ((info->t == PROJECTIVE_SUBSET) && (info->s.Holds() == VEC_SPACE)) {
  462.     return (info->d + 1 == (info->s).Dim());
  463.   } else if ((info->t == AFFINE_SUBSET) && (info->s.Holds() == PROJ_SPACE)) {
  464.     return (FALSE);
  465.   } else {
  466.     return (info->d == (info->s).Dim());
  467.   }
  468. }
  469.  
  470. // ***********************************************************************
  471. //
  472. // Returns true iff the subset is projective and was constructed
  473. // by specifying a set of base points.
  474.  
  475. Boolean SubSet::HasRemovedPoints(void)
  476. {
  477.   return ((holds == PROJECTIVE_SUBSET) && (info->substart != 0));
  478. }
  479.  
  480. // ***********************************************************************
  481. //
  482. // Returns a list of points that spans the points at infinity.  Only
  483. // applicable for affine subsets defined on projective spaces.
  484.  
  485. GeObList SubSet::AtInfinity(void)
  486. {
  487.   static char* sig = "GeObList SubSet::AtInfinity(void)";
  488.  
  489.   if ((holds != AFFINE_SUBSET) || (info->s.Holds() != PROJ_SPACE)) {
  490.     errh.ErrorExit(sig, "Must be affine subset in projective space", *this);
  491.   }   
  492.  
  493. // The tofull matrix encodes the points at infinity.  Extract the
  494. // std. coords. and build the list.
  495.  
  496.   int rows = info->tofull.Rows() - 1;
  497.   GeObList retval(rows);
  498.  
  499.   for (int i = 0; i < rows; i++) {
  500.     retval[i] = GeOb(info->accepts, info->s, info->tofull[i]);
  501.   }
  502.   return (retval);
  503. }
  504.  
  505. // ***********************************************************************
  506. //
  507. // Returns the tangent subset associated with an affine subset.
  508. // If embedding space is affine, the subset is a vector subspace of
  509. // the tangent space.  If embedding space is a vector space, the subset
  510. // is a vector subset of that space.  If it is in a projective
  511. // space, it is an error to use this.  
  512. // 
  513.  
  514. VSubSet SubSet::TangentSub(void)
  515. {
  516.   static char* sig = "VSubSet SubSet::TangentSub(void)";
  517.  
  518.   if (holds != AFFINE_SUBSET) {
  519.     errh.ErrorExit(sig, "Only allowed for affine subsets", *this);
  520.   }
  521.   if ((info->s).Holds() == PROJ_SPACE) {
  522.     errh.ErrorExit(sig, "Embedding space cannot be projective", *this);
  523.   }
  524.  
  525. // Perhaps the user has already requested the tangent space, so one has
  526. // been built.  Reuse this tangent space.
  527.  
  528.   if ((info->tansub).Holds() != NULL_SUBSET) {
  529.     return (info->tansub);
  530.   }
  531.  
  532. // If the embedding space is an affine space, and the subset is the
  533. // same dimension, the tangent subset is simply the full standard
  534. // subset of the tangent space:
  535.  
  536.   if (((info->s).Holds() == AFF_SPACE) && (this->IsFullSpace())) {
  537.     info->tansub = (info->s).GetSpace(TANGENT).FullSet();
  538.     return (info->tansub);
  539.   }
  540.  
  541. // Need to build a new tangent space.
  542.  
  543.   SubSetInfo* temp = new SubSetInfo;
  544.   Boolean isvs = (info->s).Holds() == VEC_SPACE;
  545.  
  546.   char* tag = this->BuildTagName("Subset tangent to ", *this);
  547.   strncpy(temp->name, tag, MAX_NAMESIZE - 1);
  548.   temp->name[MAX_NAMESIZE - 1] = '\0';
  549.   this->DestroyTagName(tag);
  550.  
  551.   temp->t = LINEAR_SUBSET;
  552.   temp->s = (isvs) ? info->s : (info->s).GetSpace(TANGENT);
  553.   temp->d = info->d;
  554.   temp->offset = 0.0;
  555.   temp->accepts = (isvs) ? info->accepts : AFF_VECTOR;
  556.   temp->substart = 0;
  557.  
  558. // If the parent space is a vector space, the situation is simple.
  559. // The tofull matrix omits the offset row, and we expect the
  560. // component in the offset direction to be zero.  Fromfull has
  561. // one less column:
  562.  
  563.   if (isvs) {
  564.     temp->tofull = Matrix(info->d, info->s.Dim());
  565.     for (int i = 0; i < info->d; i++) {
  566.       temp->tofull[i] = info->tofull[i];
  567.     }
  568.     temp->fullbas = info->fullbas;
  569.     temp->tstmtx = info->tstmtx;
  570.     temp->zerostart = info->isOffset;
  571.     temp->isOffset = -1;
  572.     Matrix trg = ZeroMatrix(temp->s.Dim(), temp->d);
  573.       for (int j = 0; j < temp->d; j++) {
  574.         trg[j][j] = 1.0;
  575.       }
  576.     temp->fromfull = temp->tstmtx * trg;
  577.   } else { 
  578.  
  579. // If the parent space is an affine space, the situation is more involved,
  580. // since we need to build the subset in the separate tangent vector
  581. // space.  Fortunately, the basis for the subset is expressed using
  582. // tangent space vectors and an offset, so all we have to do is
  583. // omit the offset and convert the representation:
  584.  
  585.     temp->tofull = Matrix(info->d, info->s.Dim());
  586.     Matrix convert = info->tofull * Transpose(info->s.HoistTangent());
  587.     int slot = 0;
  588.     for (int i = 0; i < convert.Rows(); i++) {
  589.       if (i != info->isOffset) {
  590.         temp->tofull[slot++] = convert[i];
  591.       }
  592.     }
  593.     Boolean errflag;
  594.     temp->fullbas = this->GSO(temp->tofull, IdentityMatrix(temp->s.Dim()),
  595.                   0, temp->s.Dim(), &errflag);
  596.     if (errflag) {
  597.       errh.ErrorExit(sig, "Unexpected error", *this);
  598.     }
  599.     temp->tstmtx = Inverse(temp->fullbas);
  600.     Matrix trg = ZeroMatrix(temp->s.Dim(), temp->d);
  601.     for (int k = 0; k < temp->d; k++) {
  602.       trg[k][k] = 1.0;
  603.     }
  604.     temp->fromfull = temp->tstmtx * trg; 
  605.     temp->zerostart = info->isOffset;
  606.     temp->isOffset = -1;
  607.   }
  608.  
  609. // Construct the new subset and return it;
  610.  
  611.   info->tansub = SubSet(temp);    
  612.   return (info->tansub);
  613. }
  614.  
  615. // ***********************************************************************
  616. //
  617. // Returns true if the given geometric object is in the subset.
  618. // (Actually, it is true if the object is within EPSILON).
  619.  
  620. Boolean SubSet::IsIn(GeOb &g)
  621. {
  622.   Scalar testval = EPSILON;
  623.  
  624. // First check if object is in the embedding space.  If not, return
  625. // false:
  626.  
  627.   GeOb temp = g.MapTo(info->accepts);
  628.   if (temp.SpaceOf() != info->s) {
  629.     return (FALSE);
  630.   }
  631.  
  632. // Then multiply the object matrix by the testing matrix, and
  633. // compare the coordinate results with the targets.
  634.  
  635.   Matrix result = temp.MatrixRep() * info->tstmtx;
  636.  
  637. // Affine subsets need a coord. == offset.  If this is an affine
  638. // subset of a projective space, we just need to make sure the
  639. // offset is not zero (but also scale the epsilon to take into account
  640. // the scaling of the offset to the desired value).  If this is a
  641. // projective subset, we need to make sure the GeOb is not a base point:
  642.  
  643.   if (info->t == AFFINE_SUBSET) {
  644.     if (info->s.Holds() == PROJ_SPACE) {
  645.       if (fabs(result[0][info->isOffset]) <= EPSILON) {
  646.         return (FALSE);
  647.       } else {
  648.         Scalar factor = fabs(info->offset / result[0][info->isOffset]);
  649.         testval = testval / factor;
  650.       }
  651.     } else {
  652.       if (fabs(result[0][info->isOffset] - info->offset) > EPSILON) {
  653.         return (FALSE);
  654.       }
  655.     }
  656.   } else if (info->t == PROJECTIVE_SUBSET) {
  657.     Boolean allzero = TRUE;
  658.     for (int i = info->substart; i < info->zerostart; i++) {
  659.       if (fabs(result[0][i]) > EPSILON) {
  660.         allzero = FALSE;
  661.     break;
  662.       }
  663.     }
  664.     if (allzero) {
  665.       return (FALSE);
  666.     }
  667.   }
  668.  
  669. // Coords. for space orthogonal to subset == 0:
  670.  
  671.   for (int i = info->zerostart; i < result.Columns(); i++) {
  672.     if (fabs(result[0][i]) > testval) {
  673.       return (FALSE);
  674.     }
  675.   }
  676.  
  677.   return (TRUE);
  678. }
  679.  
  680. // ***********************************************************************
  681. //
  682. // Given another subset, this function returns TRUE iff the given subset
  683. // is a subset of this subset.  Currently not implemented for
  684. // projective subsets with removed base points.  Note that only
  685. // comparisons of subsets of identical type are allowed.  (The system
  686. // does not test if an affine subset is a subset of a vector subset,
  687. // but simply answers FALSE).
  688.  
  689. Boolean SubSet::IsSubset(SubSet &s)
  690. {
  691.   static char* sig = "Boolean SubSet::IsSubset(SubSet&)";
  692.   int toprow;
  693.  
  694. // A quick kill is if the subsets are the same definition:
  695.  
  696.   if (info == s.info) {
  697.     return (TRUE);
  698.   }
  699.  
  700. // Check that the operation is valid:
  701.  
  702.   if (this->HasRemovedPoints() || s.HasRemovedPoints()) {
  703.     errh.ErrorExit(sig,
  704.       "Comparing projective subsets with removed points not implemented",
  705.       *this, s);
  706.   }   
  707.  
  708. // Make sure the embedding spaces and subset types match:
  709.  
  710.   if (s.info->s != info->s) {
  711.     return (FALSE);
  712.   }
  713.   if (s.holds != holds) {
  714.     return (FALSE);
  715.   }
  716.  
  717. // Hit the representation with the testing matrix:
  718.  
  719.   Matrix rep = s.ToFull() * info->tstmtx;
  720.  
  721. // For affine subsets, the origin of the argument subset frame
  722. // must be in this subset, and the tangent space elements cannot
  723. // have an offset: 
  724.  
  725.   if (s.Holds() == AFFINE_SUBSET) {
  726.     Scalar testval = EPSILON;
  727.     toprow = rep.Rows() - 1;
  728.     if (info->s.Holds() == PROJ_SPACE) {
  729.       if (fabs(rep[s.info->isOffset][info->isOffset]
  730.       * s.info->offset) <= EPSILON) {
  731.         return (FALSE);
  732.       } else {  // Scale test value to reflect normalization to offset
  733.         Scalar factor = fabs(info->offset /
  734.                  rep[s.info->isOffset][info->isOffset]);
  735.         testval = testval / factor;
  736.       }
  737.     } else {
  738.       if (fabs((rep[s.info->isOffset][info->isOffset]
  739.       * s.info->offset) - info->offset) > EPSILON) {
  740.       return (FALSE);
  741.       }
  742.     }
  743.     for (int j = info->zerostart; j < rep.Columns(); j++) {
  744.       if (fabs(rep[s.info->isOffset][j]) > testval) {
  745.         return (FALSE);
  746.       }
  747.     }
  748.     for (int i = 0; i < toprow; i++) {
  749.       if (fabs(rep[i][info->isOffset]) > EPSILON) {
  750.         return (FALSE);
  751.       }
  752.     }
  753.   } else {
  754.     toprow = rep.Rows();
  755.   }
  756.  
  757. //  Each basis vector for the argument subset or subset tangent
  758. //  space must be contained in this subset or subset tangent space:
  759.  
  760.   for (int i = 0; i < toprow; i++) {
  761.     for (int j = info->zerostart; j < rep.Columns(); j++) {
  762.       if (fabs(rep[i][j]) > EPSILON) {
  763.         return (FALSE);
  764.       }
  765.     }
  766.   }
  767.   return (TRUE);
  768. }
  769.  
  770. // ***********************************************************************
  771. // ***********************************************************************
  772. //
  773. // VSubSet Class
  774. //
  775. // ***********************************************************************
  776. // ***********************************************************************
  777. //
  778. // Private/protected member functions
  779. //
  780. // ***********************************************************************
  781. //
  782. // Builds a standard subset for a vector space.  CAUTION:  Although this
  783. // is a single-argument constructor, it is NOT intended to be used to
  784. // cast a space into a subset!
  785.  
  786. VSubSet::VSubSet(VSpace& s)
  787. {
  788. // Create the new SubSetInfo:
  789.  
  790.   info = new SubSetInfo;
  791.   type = LINEAR_SUBSET;
  792.   holds = LINEAR_SUBSET;
  793.  
  794.   char* tag = this->BuildTagName("Standard subset for ", s);
  795.   this->CopyDebugName(tag);
  796.   this->DestroyTagName(tag);
  797.  
  798.   info->s = s;
  799.   info->t = holds;
  800.   info->d = s.Dim();
  801.   info->tstmtx = IdentityMatrix(info->d);
  802.   info->fullbas = IdentityMatrix(info->d);
  803.   info->isOffset = -1;
  804.   info->substart = 0;
  805.   info->zerostart = info->d;
  806.   info->offset = 0.0;
  807.   info->fromfull = IdentityMatrix(info->d);
  808.   info->tofull = IdentityMatrix(info->d);
  809.   info->accepts = s.NativeType();
  810. }
  811.  
  812. // ***********************************************************************
  813. //
  814. // Public member functions
  815. //
  816. // ***********************************************************************
  817. //
  818. // Memberwise assignment:
  819. //
  820.  
  821. VSubSet& VSubSet::operator=(VSubSet &s)
  822. {
  823.   info = s.info;
  824.   holds = s.holds;
  825.   return *this;
  826. }
  827.  
  828. // ***********************************************************************
  829. //
  830. // A subset can only be cast down to a vector subset if it is one (no
  831. // conversions).
  832. //
  833.  
  834. VSubSet::VSubSet(SubSet& s) : (s)
  835. {
  836.   type = LINEAR_SUBSET;
  837.   if ((holds != LINEAR_SUBSET) && (holds != NULL_SUBSET)) {
  838.       errh.ErrorExit("VSubSet::VSubSet(SubSet&)",
  839.                 "Attempted to cast a non-vector subset to a vector subset", s);
  840.   }
  841. }
  842.  
  843. // ***********************************************************************
  844. //
  845. // Output for debugging:
  846. //
  847.  
  848. void VSubSet::debug_out(ostream &c, int indent) 
  849. {
  850.   static char* name = "VSubSet object\n";
  851.   this->data_out(c, indent, name);
  852.   return;
  853. }
  854.  
  855. // ***********************************************************************
  856. //
  857. //  Build a vector subset of a vector space.
  858. //
  859.  
  860. VSubSet::VSubSet(char* namein, VSpace &s, GeObList &v)
  861. {
  862.   static char* sig = "VSubSet::VSubSet(char*, VSpace&, GeObList&)";
  863.   Boolean errflag;
  864.  
  865.   if (v.Length() < 1) {
  866.     errh.ErrorExit(sig, "Incorrect number of points specified", v);
  867.   }
  868.  
  869. // Create the SubSetInfo block:
  870.  
  871.   info = new SubSetInfo;
  872.   type = LINEAR_SUBSET;
  873.   holds = LINEAR_SUBSET;
  874.  
  875. // Important info:
  876.  
  877.   info->t = holds; 
  878.   info->d = v.Length();
  879.   info->s = s;
  880.   info->accepts = s.NativeType();
  881.   info->isOffset = -1;
  882.   info->substart = 0;
  883.   info->zerostart = info->d;
  884.   info->offset = 0.0;
  885.   this->CopyDebugName(namein);
  886.  
  887. // Make sure all the spanning objects are in the specified space:
  888.  
  889.   Matrix spantemp(v.Length(), s.Dim());
  890.   if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
  891.     errh.ErrorExit(sig,
  892.            "Object cannot be mapped into specified space", v, s);
  893.   }
  894.  
  895. // The user provides a list of geometric objects that span
  896. // some n-dimensional portion of the m-dimensional space.  They must
  897. // be independent.  We need to construct an orthonormal basis for
  898. // the embedding space such that the first n vectors span the 
  899. // subset, and the remaining vectors are orthogonal to the subset.
  900. // The orthonormal basis is built using Gram-Schmidt othogonalization.
  901. // Normalize the first vector:
  902.  
  903.   Scalar mag = sqrt((spantemp[0] * Transpose(spantemp[0]))[0][0]);
  904.   Matrix have = spantemp[0] / mag;
  905.  
  906. // Gram-Schmidt:
  907.  
  908.   info->tofull = this->GSO(have, spantemp, 1, info->d, &errflag);
  909.   if (errflag) {
  910.     errh.ErrorExit(sig, "Objects to define subset must be independent", v, s);
  911.   }
  912.   info->fullbas = this->GSO(info->tofull, IdentityMatrix(s.Dim()),
  913.                 0, s.Dim(), &errflag);
  914.   if (errflag) {
  915.     errh.ErrorExit(sig, "Unexpected error", v, s);
  916.   }
  917.  
  918. // Build the testing matrix and the projection matrix to convert 
  919. // "close" full space representations to subset representations:
  920.  
  921.   info->tstmtx = Inverse(info->fullbas);
  922.   Matrix trg = ZeroMatrix(s.Dim(), info->d);
  923.   for (int i = 0; i < info->d; i++) {
  924.     trg[i][i] = 1.0;
  925.   }
  926.   info->fromfull = info->tstmtx * trg; 
  927. } 
  928.  
  929. // ***********************************************************************
  930. // ***********************************************************************
  931. //
  932. // ASubSet Class
  933. //
  934. // ***********************************************************************
  935. // ***********************************************************************
  936. //
  937. // Private/protected member functions
  938. //
  939. // ***********************************************************************
  940. //
  941. // Builds a standard affine subset for a space.  CAUTION:  Although this
  942. // is a single-argument constructor, it is NOT intended to be used to
  943. // cast a space into a subset!
  944. //
  945.  
  946. ASubSet::ASubSet(Space& s)
  947. {
  948.   int size = s.MatrixSize();
  949.  
  950. // Create the new SubSetInfo:
  951.  
  952.   info = new SubSetInfo;
  953.   type = AFFINE_SUBSET;
  954.   holds = AFFINE_SUBSET;
  955.  
  956.   char* tag = this->BuildTagName("Standard affine subset for ", s);
  957.   this->CopyDebugName(tag);
  958.   this->DestroyTagName(tag);
  959.  
  960.   info->t = holds;
  961.   info->d = (s.Holds() == VEC_SPACE) ? s.Dim() - 1 : s.Dim();
  962.   info->s = s;
  963.   info->tstmtx = IdentityMatrix(size);
  964.   info->fullbas = IdentityMatrix(size);
  965.   info->isOffset = size - 1;
  966.   info->substart = 0;
  967.   info->zerostart = size;
  968.   info->offset = 1.0;
  969.   info->fromfull = IdentityMatrix(size);
  970.   info->tofull = IdentityMatrix(size);
  971.   info->accepts = s.NativeType();
  972. }
  973.  
  974. // ***********************************************************************
  975. //
  976. // Used to create a consistent standard affine subset, given a
  977. // projective map.
  978. //
  979.  
  980. ASubSet::ASubSet(ASubSet& a, ProjectiveMap& m)
  981. {
  982.   static char* sig = "ASubSet::ASubSet(ASubSet&, ProjectiveMap&)";
  983.   Boolean errflag;
  984.  
  985. // Create the new SubSetInfo block:
  986.  
  987.   info = new SubSetInfo;
  988.   type = AFFINE_SUBSET;
  989.   holds = AFFINE_SUBSET;
  990.  
  991.   info->t = holds;
  992.   info->s = m.Range().EmbeddingSpace();
  993.  
  994.   int dim = info->s.MatrixSize();
  995.  
  996.   info->d = (info->s.Holds() == VEC_SPACE) ? info->s.Dim() - 1 : info->s.Dim();
  997.   info->accepts = info->s.NativeType();
  998.   info->isOffset = dim - 1;
  999.   info->substart = 0;
  1000.   info->zerostart = dim;
  1001.  
  1002.   char* tag = this->BuildTagName("Standard affine subset for ", info->s);
  1003.   this->CopyDebugName(tag);
  1004.   this->DestroyTagName(tag);
  1005.  
  1006. // The map had better be invertible.  Send the subset basis through
  1007. // the map, and recreate the subset:
  1008.  
  1009.   if (!m.Invertible()) {
  1010.     errh.ErrorExit(sig, "Requires invertible map", a, m);
  1011.   }
  1012.  
  1013.   if (a.EmbeddingSpace() != m.Domain().EmbeddingSpace()) {
  1014.     errh.ErrorExit(sig, "Domain mismatch", a, m);
  1015.   }
  1016.  
  1017. // This approach works, but is not elegant.
  1018. // Map the offset vector as a linear functional to keep it
  1019. // perpendicular to tangent space.  Map the subset basis
  1020. // as usual:
  1021.  
  1022.   Matrix pt = ((a.AugBasis())[a.IsOffset()] *
  1023.           Inverse(Transpose(m.MatrixRep())));
  1024.   Scalar mag = sqrt((pt * Transpose(pt))[0][0]);
  1025.   info->offset = a.Offset() / mag;
  1026.  
  1027.   Matrix span2 = a.AugBasis() * m.MatrixRep();
  1028.  
  1029. // Keep the offset vector intact by making it the first vector
  1030. // in the new basis:
  1031.  
  1032.   span2[a.IsOffset()] = span2[0];
  1033.   span2[0] = pt[0];
  1034.  
  1035. // Build a set of orthonormal vectors for the subspace, then
  1036. // augment with vectors orthogonal to the subspace:
  1037.  
  1038.   Matrix have = span2[0] / mag;
  1039.   info->tofull = this->GSO(have, span2, 1, dim, &errflag);
  1040.   if (errflag) {
  1041.     errh.ErrorExit(sig, "Unexpected error A", a, m);
  1042.   }
  1043.  
  1044. // Swap unchanged offset vector back into position:
  1045.  
  1046.   Matrix hold = info->tofull[0];
  1047.   info->tofull[0] = info->tofull[info->isOffset];
  1048.   info->tofull[info->isOffset] = hold[0];
  1049.  
  1050.   info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
  1051.                 0, dim, &errflag);
  1052.   if (errflag) {
  1053.     errh.ErrorExit(sig, "Unexpected error B", a, m);
  1054.   }
  1055.  
  1056. // Calculate the testing matrix:
  1057.  
  1058.   info->tstmtx = Inverse(info->fullbas);
  1059.  
  1060. // Calculate fromfull:
  1061.  
  1062.   Matrix trg = ZeroMatrix(dim, dim);
  1063.   for (int i = 0; i < dim; i++) {
  1064.     trg[i][i] = 1.0;
  1065.   }
  1066.   info->fromfull = info->tstmtx * trg; 
  1067. }
  1068.  
  1069. // ***********************************************************************
  1070. //
  1071. // Public member functions
  1072. //
  1073. // ***********************************************************************
  1074. //
  1075. // A subset can only be cast down to an affine subset if it is one (no
  1076. // conversions).
  1077.  
  1078. ASubSet::ASubSet(SubSet& s) : (s)
  1079. {
  1080.   type = AFFINE_SUBSET;
  1081.   if ((holds != AFFINE_SUBSET) && (holds != NULL_SUBSET)) {
  1082.       errh.ErrorExit("ASubSet::ASubSet(SubSet&)",
  1083.                "Attempted to cast a non-affine subset to an affine subset", s);
  1084.   }
  1085. }
  1086.  
  1087. // ***********************************************************************
  1088. //
  1089. // Memberwise assignment:
  1090. //
  1091.  
  1092. ASubSet& ASubSet::operator=(ASubSet& s)
  1093. {
  1094.   info = s.info;
  1095.   holds = s.holds;
  1096.   return *this;
  1097. }  
  1098.  
  1099. // ***********************************************************************
  1100. //
  1101. // Output for debugging:
  1102. //
  1103.  
  1104. void ASubSet::debug_out(ostream &c, int indent) 
  1105. {
  1106.   static char* name = "ASubSet object\n";
  1107.   this->data_out(c, indent, name);
  1108.   return;
  1109. }
  1110.  
  1111. // ***********************************************************************
  1112. //
  1113. //  Build an affine subset of a vector or affine space.
  1114. //
  1115.  
  1116. ASubSet::ASubSet(char* namein, Space &s, GeObList &v)
  1117. {
  1118.   static char* sig = "ASubSet::ASubSet(char*, Space&, GeObList&)";
  1119.   Boolean errflag;
  1120.  
  1121.   if ((s.Holds() != VEC_SPACE) && (s.Holds() != AFF_SPACE)) {
  1122.     errh.ErrorExit(sig, "Vector or affine space required", s);
  1123.   }
  1124.  
  1125.   if (v.Length() < 1) {
  1126.     errh.ErrorExit(sig, "Incorrect number of points specified", v);
  1127.   }
  1128.  
  1129. // Create the new SubSetInfo block:
  1130.  
  1131.   info = new SubSetInfo;
  1132.   type = AFFINE_SUBSET;
  1133.   holds = AFFINE_SUBSET;
  1134.  
  1135. // Important info:
  1136.  
  1137.   info->t = holds; 
  1138.   info->d = v.Length() - 1;
  1139.   info->s = s;
  1140.   info->accepts = s.NativeType();
  1141.   info->isOffset = info->d;
  1142.   info->substart = 0;
  1143.   info->zerostart = info->d + 1;
  1144.   this->CopyDebugName(namein);
  1145.  
  1146. // Make sure all the spanning objects are in the specified space:
  1147.  
  1148.   int dim = s.MatrixSize();
  1149.   int num = v.Length();
  1150.  
  1151.   Matrix spantemp(num, dim);
  1152.   if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
  1153.     errh.ErrorExit(sig, "Object cannot be mapped into specified space", v, s);
  1154.   }
  1155.  
  1156. // The subset basis will be a frame for the affine subset, so we need
  1157. // to derive vectors in the tangent subset.  Keep the last point
  1158. // in the list as the point origin for the frame, and keep a separate
  1159. // copy of the point's matrix rep. to use for finding the offset.
  1160.  
  1161.   Matrix span2(num, dim);
  1162.   Matrix pt = spantemp[num - 1];
  1163.   int slot = 0;
  1164.   for (int i = 0; i < num; i++) {
  1165.     if (i == num - 1) {
  1166.       span2[i] = spantemp[i];
  1167.     } else {
  1168.  
  1169. // The following line of code produces a bug with g++ 1.37.1
  1170. // (no problem with CFRONT 1.2) when using dynamically allocated matrices:
  1171. //    span2[slot++] = (spantemp[i] - pt)[0];
  1172. // Replace it with:
  1173.  
  1174.       Matrix temp = (spantemp[i] - pt[0]);
  1175.       span2[slot++] = temp[0];
  1176.     }
  1177.   }
  1178.  
  1179. // Build a set of orthonormal vectors for the subspace, then
  1180. // augment with vectors orthogonal to the subspace:
  1181.  
  1182.   Scalar mag = sqrt((span2[0] * Transpose(span2[0]))[0][0]);
  1183.   Matrix have = span2[0] / mag;
  1184.   info->tofull = this->GSO(have, span2, 1, num, &errflag);
  1185.   if (errflag) {
  1186.     errh.ErrorExit(sig, "Objects to define subset must be independent", v, s);
  1187.   }
  1188.  
  1189.   info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
  1190.                 0, dim, &errflag);
  1191.   if (errflag) {
  1192.     errh.ErrorExit(sig, "Unexpected error", v, s);
  1193.   }
  1194.  
  1195. // Testing matrix:
  1196.  
  1197.   info->tstmtx = Inverse(info->fullbas);
  1198.  
  1199. // Record the offset of the affine hyperplane (the projection of
  1200. // the original frame origin point onto the offset dimension):
  1201.  
  1202.   info->offset = (pt * Transpose(info->tofull[info->isOffset]))[0][0];
  1203.  
  1204. // Calculate fromfull:
  1205.  
  1206.   Matrix trg = ZeroMatrix(dim, num);
  1207.   for (int j = 0; j < num; j++) {
  1208.     trg[j][j] = 1.0;
  1209.   }
  1210.   info->fromfull = info->tstmtx * trg; 
  1211. } 
  1212.  
  1213. // ***********************************************************************
  1214. //
  1215. // Build an affine subset of a projective space. The key here is 
  1216. // that the combination of the "regular" point and the
  1217. // "points at infinity" specify a projective subspace;
  1218. // subtracting the points at infinity turns this into an
  1219. // affine subset.
  1220. //
  1221.  
  1222. ASubSet::ASubSet(char* namein, PSpace& s, GeOb& v, GeObList& inf)
  1223. {
  1224.   static char* sig = "ASubSet::ASubSet(char*, PSpace&, GeOb&, GeObList&)";
  1225.   Boolean errflag;
  1226.  
  1227. // The number of points at infinity that are provided needs to be
  1228. // appropriate.  For an m-dimensional affine subset, we need m+1
  1229. // projective points, of which m are designated to be points at
  1230. // infinity, and only 1 is a "regular" point.
  1231.  
  1232.   if (inf.Length() < 1) {
  1233.     errh.ErrorExit(sig, "Incorrect number of points specified", inf);
  1234.   }
  1235.  
  1236. // Create a new SubSetInfo block:
  1237.  
  1238.   info = new SubSetInfo;
  1239.   type = AFFINE_SUBSET;
  1240.   holds = AFFINE_SUBSET;
  1241.  
  1242. // Important info:
  1243.  
  1244.   info->t = holds; 
  1245.   info->d = inf.Length();
  1246.   info->s = s;
  1247.   info->accepts = s.NativeType();
  1248.   info->isOffset = info->d;
  1249.   info->substart = 0;
  1250.   info->zerostart = info->d + 1;
  1251.   this->CopyDebugName(namein);
  1252.  
  1253. // Make sure all the spanning objects are in the specified space:
  1254.  
  1255.   int dim = s.MatrixSize();
  1256.  
  1257.   Matrix vtemp(1, dim);
  1258.   Matrix inftemp(inf.Length(), dim);
  1259.   GeObList vlst(v);
  1260.  
  1261.   if (!(this->ConvertList(vlst, info->accepts, s, &vtemp))) {
  1262.     errh.ErrorExit(sig, "Object cannot be mapped into specified space", v, s);
  1263.   }
  1264.   if (!(this->ConvertList(inf, info->accepts, s, &inftemp))) {
  1265.     errh.ErrorExit(sig,
  1266.            "Object cannot be mapped into specified space", inf, s);
  1267.   }
  1268.  
  1269. // Combine the matrix reps. of the two lists into a single matrix.
  1270.  
  1271.   int num = info->d + 1;
  1272.   Matrix span(num, dim);
  1273.   span[num - 1] = vtemp[0];
  1274.   for (int i = 0; i < num - 1; i++) {
  1275.     span[i] = inftemp[i];
  1276.   }
  1277.  
  1278. // Build a set of orthonormal vectors for the subset, then
  1279. // augment with vectors orthogonal to the subset:
  1280.  
  1281.   Scalar mag = sqrt((span[0] * Transpose(span[0]))[0][0]);
  1282.   Matrix have = span[0] / mag;
  1283.   info->tofull = this->GSO(have, span, 1, num, &errflag);
  1284.   if (errflag) {
  1285.     errh.ErrorExit(sig,
  1286.            "Objects to define subset must be independent", v, inf);
  1287.   }
  1288.   info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
  1289.                 0, dim, &errflag);
  1290.   if (errflag) {
  1291.     errh.ErrorExit(sig, "Unexpected error", v, inf);
  1292.   }
  1293.  
  1294. // Build the testing matrix:
  1295.  
  1296.   info->tstmtx = Inverse(info->fullbas);
  1297.  
  1298. // Set the offset of the affine hyperplane at 1.0 (it can be any
  1299. // arbitrary non-zero value for affine subsets of projective spaces):
  1300.  
  1301.   info->offset = 1.0;
  1302.  
  1303. // Calculate fromfull:
  1304.  
  1305.   Matrix trg = ZeroMatrix(dim, num);
  1306.   for (int j = 0; j < num; j++) {
  1307.     trg[j][j] = 1.0;
  1308.   }
  1309.   info->fromfull = info->tstmtx * trg; 
  1310. } 
  1311.  
  1312. // ***********************************************************************
  1313. // ***********************************************************************
  1314. //
  1315. // PSubSet Class
  1316. //
  1317. // ***********************************************************************
  1318. // ***********************************************************************
  1319. //
  1320. // Private/protected member functions
  1321. //
  1322. // ***********************************************************************
  1323. //
  1324. // Builds a standard projective subset for a space.
  1325. //
  1326.  
  1327. PSubSet::PSubSet(Space& s)
  1328. {
  1329.   static char* sig = "PSubSet::PSubSet(Space&)";
  1330.   int size = s.MatrixSize();
  1331.  
  1332.   if ((s.Holds() != VEC_SPACE) && (s.Holds() != PROJ_SPACE)) {
  1333.     errh.ErrorExit(sig, "Vector or projective space required", s);
  1334.   }
  1335.  
  1336. // Create a new SubSetInfo block:
  1337.  
  1338.   info = new SubSetInfo;
  1339.   type = PROJECTIVE_SUBSET;
  1340.   holds = PROJECTIVE_SUBSET;
  1341.  
  1342.   char* tag = this->BuildTagName("Standard proj. subset for ", s);
  1343.   this->CopyDebugName(tag);
  1344.   this->DestroyTagName(tag);
  1345.  
  1346.   info->t = holds;
  1347.   info->d = (s.Holds() == VEC_SPACE) ? s.Dim() - 1 : s.Dim();
  1348.   info->s = s;
  1349.   info->tstmtx = IdentityMatrix(size);
  1350.   info->fullbas = IdentityMatrix(size);
  1351.   info->isOffset = -1;
  1352.   info->substart = 0;
  1353.   info->zerostart = size;
  1354.   info->offset = 0.0;
  1355.   info->fromfull = IdentityMatrix(size);
  1356.   info->tofull = IdentityMatrix(size);
  1357.  
  1358.   info->accepts = s.NativeType();
  1359.   if (info->accepts == VECTOR) {
  1360.     info->accepts = VEC_EC;
  1361.   } else if (info->accepts == AFF_VECTOR) {
  1362.     info->accepts = AFF_VEC_EC;
  1363.   }
  1364. }
  1365.  
  1366. // ***********************************************************************
  1367. //
  1368. // Public member functions
  1369. //
  1370. // ***********************************************************************
  1371. //
  1372. // A subset can only be cast down to a proj. subset if it is one (no
  1373. // conversions).
  1374.  
  1375. PSubSet::PSubSet(SubSet& s) : (s)
  1376. {
  1377.   type = PROJECTIVE_SUBSET;
  1378.   if ((holds != PROJECTIVE_SUBSET) && (holds != NULL_SUBSET)) {
  1379.       errh.ErrorExit("PSubSet::PSubSet(SubSet&)",
  1380.        "Attempted to cast a non-projective subset to a projective subset", s);
  1381.   }
  1382. }
  1383.  
  1384. // ***********************************************************************
  1385. //
  1386. // Memberwise assignment:
  1387. //
  1388.  
  1389. PSubSet& PSubSet::operator=(PSubSet &s)
  1390. {
  1391.   info = s.info;
  1392.   holds = s.holds;
  1393.   return *this;
  1394. }  
  1395.  
  1396. // ***********************************************************************
  1397.  
  1398. void PSubSet::debug_out(ostream &c, int indent) 
  1399. {
  1400.   static char* name = "PSubSet object\n";
  1401.   this->data_out(c, indent, name);
  1402.   return;
  1403. }
  1404.  
  1405. // ***********************************************************************
  1406. //
  1407. //  Build a projective subspace, without base point specification.
  1408. //
  1409.  
  1410. PSubSet::PSubSet(char* namein, Space &s, GeObList &v)
  1411. {
  1412.   static char* sig = "PSubSet::PSubSet(char*, Space&, GeObList&)";
  1413.   Boolean errflag;
  1414.  
  1415.   if ((s.Holds() != VEC_SPACE) && (s.Holds() != PROJ_SPACE)) {
  1416.     errh.ErrorExit(sig, "Vector or Projective space required", s);
  1417.   }
  1418.  
  1419.   if (v.Length() < 1) {
  1420.     errh.ErrorExit(sig, "Incorrect number of points specified", v);
  1421.   }
  1422.  
  1423. // Create a new SubSetInfo block:
  1424.  
  1425.   info = new SubSetInfo;
  1426.   type = PROJECTIVE_SUBSET;
  1427.   holds = PROJECTIVE_SUBSET;
  1428.  
  1429. // Important info:
  1430.  
  1431.   info->t = holds; 
  1432.   info->d = v.Length() - 1;
  1433.   info->s = s;
  1434.   info->isOffset = -1;
  1435.   info->substart = 0;
  1436.   info->zerostart = info->d + 1;
  1437.   info->offset = 0.0;
  1438.   this->CopyDebugName(namein);
  1439.  
  1440.   info->accepts = s.NativeType();
  1441.   if (info->accepts == VECTOR) {
  1442.     info->accepts = VEC_EC;
  1443.   } else if (info->accepts == AFF_VECTOR) {
  1444.     info->accepts = AFF_VEC_EC;
  1445.   }
  1446.  
  1447. // Make sure all the spanning objects are in the specified space:
  1448.  
  1449.   int dim = s.MatrixSize();
  1450.   int num = v.Length();
  1451.  
  1452.   Matrix spantemp(num, dim);
  1453.   if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
  1454.     errh.ErrorExit(sig,
  1455.            "Object cannot be mapped into specified space", v, s);
  1456.   }
  1457.  
  1458. // Build a set of orthonormal vectors for the subspace, then
  1459. // augment with vectors orthogonal to the subspace:
  1460.  
  1461.   Scalar mag = sqrt((spantemp[0] * Transpose(spantemp[0]))[0][0]);
  1462.   Matrix have = spantemp[0] / mag;
  1463.   info->tofull = this->GSO(have, spantemp, 1, num, &errflag);
  1464.   if (errflag) {
  1465.     errh.ErrorExit(sig, "Objects to define subset must be independent", v, s);
  1466.   }
  1467.   info->fullbas = this->GSO(info->tofull, IdentityMatrix(dim),
  1468.                 0, dim, &errflag);
  1469.   if (errflag) {
  1470.     errh.ErrorExit(sig, "Unexpected error", v, s);
  1471.   }
  1472.  
  1473. // Build testing matrix and fromfull matrix:
  1474.  
  1475.   info->tstmtx = Inverse(info->fullbas);
  1476.   Matrix trg = ZeroMatrix(dim, num);
  1477.   for (int i = 0; i < num; i++) {
  1478.     trg[i][i] = 1.0;
  1479.   }
  1480.   info->fromfull = info->tstmtx * trg; 
  1481. } 
  1482.  
  1483. // ***********************************************************************
  1484. //
  1485. // Build a projective subset of a projective space by selection of
  1486. // base points.  In combination, base points and other points
  1487. // define a projective subset.  The subsequent removal of the
  1488. // base points reduces the dimensionality of the subset.  The
  1489. // resulting dimensionality equals that of the corresponding
  1490. // primitive geometric form (e.g. pencil or bundle).
  1491. //
  1492.  
  1493. PSubSet::PSubSet(char* namein, Space& s, GeObList& bp, GeObList& v)
  1494. {
  1495.   static char* sig = "PSubSet::PSubSet(char*, Space&, GeObList&, GeObList&)";
  1496.   Boolean errflag;
  1497.   int num = v.Length() + bp.Length();
  1498.  
  1499.   if ((s.Holds() != VEC_SPACE) && (s.Holds() != PROJ_SPACE)) {
  1500.     errh.ErrorExit(sig, "Vector or Projective space required", s);
  1501.   }
  1502.  
  1503.   if ((v.Length() < 1) || (bp.Length() < 1)) {
  1504.     errh.ErrorExit(sig, "Incorrect number of points specified", bp, v);
  1505.   }
  1506.  
  1507. // Create a new SubSetInfo block:
  1508.  
  1509.   info = new SubSetInfo;
  1510.   type = PROJECTIVE_SUBSET;
  1511.   holds = PROJECTIVE_SUBSET;
  1512.  
  1513. // Important info:
  1514.  
  1515.   info->t = holds; 
  1516.   info->d = v.Length() - 1;
  1517.   info->s = s;
  1518.   info->isOffset = -1;
  1519.   info->substart = bp.Length();
  1520.   info->zerostart = num;
  1521.   info->offset = 0.0;
  1522.   this->CopyDebugName(namein);
  1523.  
  1524.   info->accepts = s.NativeType();
  1525.   if (info->accepts == VECTOR) {
  1526.     info->accepts = VEC_EC;
  1527.   } else if (info->accepts == AFF_VECTOR) {
  1528.     info->accepts = AFF_VEC_EC;
  1529.   }
  1530.  
  1531. // Make sure all the spanning objects are in the specified space:
  1532.  
  1533.   int dim = s.MatrixSize();
  1534.  
  1535.   Matrix spantemp(v.Length(), dim);
  1536.   Matrix bptemp(bp.Length(), dim);
  1537.  
  1538.   if (!(this->ConvertList(v, info->accepts, s, &spantemp))) {
  1539.     errh.ErrorExit(sig,
  1540.            "Object cannot be mapped into specified space", v, s);
  1541.   }
  1542.   if (!(this->ConvertList(bp, info->accepts, s, &bptemp))) {
  1543.     errh.ErrorExit(sig,
  1544.            "Object cannot be mapped into specified space", bp, s);
  1545.   }
  1546.  
  1547. // Build a set of orthonormal vectors that span the space of base points.
  1548. // Then build an orthogonal basis to represent the subset, and finally
  1549. // add enough vectors to fill the basis out.  Note that we 
  1550. // CANNOT build a matrix that converts basis std. coords. to full space
  1551. // std. coords (tofull), since fromfull is not invertible.
  1552.  
  1553.   Scalar mag = sqrt((bptemp[0] * Transpose(bptemp[0]))[0][0]);
  1554.   Matrix have = bptemp[0] / mag;
  1555.   Matrix temp = this->GSO(have, bptemp, 1, bp.Length(), &errflag);
  1556.   if (errflag) {
  1557.     errh.ErrorExit(sig,
  1558.            "Objects to define subset must be independent", s, bp, v);
  1559.   }
  1560.   temp = this->GSO(temp, spantemp, 0, num, &errflag);
  1561.   if (errflag) {
  1562.     errh.ErrorExit(sig,
  1563.            "Objects to define subset must be independent", s, bp, v);
  1564.   }
  1565.   info->fullbas = this->GSO(temp, IdentityMatrix(dim), 0, dim, &errflag);
  1566.   if (errflag) {
  1567.     errh.ErrorExit(sig, "Unexpected error", s, bp, v);
  1568.   }
  1569.  
  1570. // Build testing matrix and fromfull matrix:
  1571.  
  1572.   info->tstmtx = Inverse(info->fullbas);
  1573.  
  1574.   Matrix trg = ZeroMatrix(dim, v.Length());
  1575.   int slot = 0;
  1576.   for (int i = bp.Length(); i < num; i++) {
  1577.     trg[i][slot++] = 1.0;
  1578.   }
  1579.   info->fromfull = info->tstmtx * trg; 
  1580. } 
  1581.  
  1582.  
  1583. // ***********************************************************************
  1584. // ***********************************************************************
  1585. //
  1586. // Create the error information block:
  1587. //
  1588.  
  1589. static SubSetInfo errInfo("Uninitialized subset", 0);
  1590.  
  1591. // ***********************************************************************
  1592.